/*
Copyright (c) 2009-2010, Warwick Warp Limited
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of the Warwick Warp Limited nor the
      names of its contributors may be used to endorse or promote products
      derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL WARWICK WARP LIMITED BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#include "warray_blas.h"

#if defined(USE_ACCELERATE) || defined(USE_MKL) || defined(USE_CBLAS)

#if defined(USE_MKL)
#include "pmkl.h"
#elif defined(USE_ACCELERATE)
#include <Accelerate/Accelerate.h>
#elif defined(USE_CBLAS)
#include "cblas.h"
#endif

namespace {

inline CBLAS_TRANSPOSE to_blas_trans(warp::blas::conjtrans_e v)
{
	switch(v)
	{
	default:
	case warp::blas::dont_transpose:
		return CblasNoTrans;
	case warp::blas::transpose:
		return CblasTrans;
	case warp::blas::conjugate_transpose:
		return CblasConjTrans;
	}
}

inline CBLAS_UPLO to_blas_uplo(warp::blas::uplo_e v)
{
	switch(v)
	{
	default:
	case warp::blas::upper_triangular:
		return CblasUpper;
	case warp::blas::lower_triangular:
		return CblasLower;
	}
}


} // namespace

template<> float warp::blas::asum<float>(size_t n, const float *x, int incx)
{
	return cblas_sasum(n,x,incx);
}

template<> double warp::blas::asum<double>(size_t n, const double *x, int incx)
{
	return cblas_dasum(n,x,incx);
}

template<> float warp::blas::asum<float>(size_t n, const std::complex<float> *x, int incx)
{
	return cblas_scasum(n,x,incx);
}

template<> double warp::blas::asum<double>(size_t n, const std::complex<double> *x, int incx)
{
	return cblas_dzasum(n,x,incx);
}

template<> void warp::blas::copy<float>(size_t n, const float *x, int incx, float *y, int incy)
{
	cblas_scopy(n,x,incx,y,incy);
}

template<> void warp::blas::copy<double>(size_t n, const double *x, int incx, double *y, int incy)
{
	cblas_dcopy(n,x,incx,y,incy);
}

template<> void warp::blas::copy< std::complex<float> >(size_t n, const std::complex<float> *x, int incx, std::complex<float> *y, int incy)
{
	cblas_ccopy(n,x,incx,y,incy);
}

template<> void warp::blas::copy< std::complex<double> >(size_t n, const std::complex<double> *x, int incx, std::complex<double> *y, int incy)
{
	cblas_zcopy(n,x,incx,y,incy);
}

template<> float warp::blas::dot<float>(size_t n, const float *x, int incx, const float *y, int incy)
{
	return cblas_sdot(n,x,incx,y,incy);
}

template<> double warp::blas::dot<double>(size_t n, const double *x, int incx, const double *y, int incy)
{
	return cblas_ddot(n,x,incx,y,incy);
}

template<> std::complex<float> warp::blas::dot< std::complex<float> >(size_t n, const std::complex<float> *x, int incx, const std::complex<float> *y, int incy)
{
	std::complex<float> ans;
	cblas_cdotu_sub(n,x,incx,y,incy,&ans);
	return ans;
}

template<> std::complex<double> warp::blas::dot< std::complex<double> >(size_t n, const std::complex<double> *x, int incx, const std::complex<double> *y, int incy)
{
	std::complex<double> ans;
	cblas_zdotu_sub(n,x,incx,y,incy,&ans);
	return ans;
}

template<> std::complex<float> warp::blas::dotc<float>(size_t n, const std::complex<float> *x, int incx, const std::complex<float> *y, int incy)
{
	std::complex<float> ans;
	cblas_cdotc_sub(n,x,incx,y,incy,&ans);
	return ans;
}

template<> std::complex<double> warp::blas::dotc<double>(size_t n, const std::complex<double> *x, int incx, const std::complex<double> *y, int incy)
{
	std::complex<double> ans;
	cblas_zdotc_sub(n,x,incx,y,incy,&ans);
	return ans;
}

template<> float warp::blas::nrm2<float>(size_t n, const float *x, int incx)
{
	return cblas_snrm2(n,x,incx);
}

template<> double warp::blas::nrm2<double>(size_t n, const double *x, int incx)
{
	return cblas_dnrm2(n,x,incx);
}

template<> float warp::blas::nrm2<float>(size_t n, const std::complex<float> *x, int incx)
{
	return cblas_scnrm2(n,x,incx);
}

template<> double warp::blas::nrm2<double>(size_t n, const std::complex<double> *x, int incx)
{
	return cblas_dznrm2(n,x,incx);
}

template<> void warp::blas::rot<float>(size_t n, float *x, int incx, float *y, int incy, const float& c, const float& s)
{
	cblas_srot(n,x,incx,y,incy,c,s);
}

template<> void warp::blas::rot<double>(size_t n, double *x, int incx, double *y, int incy, const double& c, const double& s)
{
	cblas_drot(n,x,incx,y,incy,c,s);
}

template<> void warp::blas::scal<float>(size_t n, const float& a, float *x, int incx)
{
	cblas_sscal(n,a,x,incx);
}

template<> void warp::blas::scal<double>(size_t n, const double& a, double *x, int incx)
{
	cblas_dscal(n,a,x,incx);
}

template<> void warp::blas::scal< std::complex<float> >(size_t n, const std::complex<float>& a, std::complex<float> *x, int incx)
{
	cblas_cscal(n,&a,x,incx);
}

template<> void warp::blas::scal< std::complex<double> >(size_t n, const std::complex<double>& a, std::complex<double> *x, int incx)
{
	cblas_zscal(n,&a,x,incx);
}

template<> void warp::blas::scal<float>(size_t n, const float& a, std::complex<float> *x, int incx)
{
	cblas_csscal(n,a,x,incx);
}

template<> void warp::blas::scal<double>(size_t n, const double& a, std::complex<double> *x, int incx)
{
	cblas_zdscal(n,a,x,incx);
}

template<> void warp::blas::swap<float>(size_t n, float *x, int incx, float *y, int incy)
{
	cblas_sswap(n,x,incx,y,incy);
}

template<> void warp::blas::swap<double>(size_t n, double *x, int incx, double *y, int incy)
{
	cblas_dswap(n,x,incx,y,incy);
}

template<> void warp::blas::swap< std::complex<float> >(size_t n, std::complex<float> *x, int incx, std::complex<float> *y, int incy)
{
	cblas_cswap(n,x,incx,y,incy);
}

template<> void warp::blas::swap< std::complex<double> >(size_t n, std::complex<double> *x, int incx, std::complex<double> *y, int incy)
{
	cblas_zswap(n,x,incx,y,incy);
}

template<> ptrdiff_t warp::blas::amax<float>(size_t n, const float *x, int incx)
{
	return cblas_isamax(n,x,incx);
}

template<> ptrdiff_t warp::blas::amax<double>(size_t n, const double *x, int incx)
{
	return cblas_idamax(n,x,incx);
}

template<> ptrdiff_t warp::blas::amax<float>(size_t n, const std::complex<float> *x, int incx)
{
	return cblas_icamax(n,x,incx);
}

template<> ptrdiff_t warp::blas::amax<double>(size_t n, const std::complex<double> *x, int incx)
{
	return cblas_izamax(n,x,incx);
}

#ifndef USE_ACCELERATE
template<> ptrdiff_t warp::blas::amin<float>(size_t n, const float *x, int incx)
{
	return cblas_isamin(n,x,incx);
}

template<> ptrdiff_t warp::blas::amin<double>(size_t n, const double *x, int incx)
{
	return cblas_idamin(n,x,incx);
}

template<> ptrdiff_t warp::blas::amin<float>(size_t n, const std::complex<float> *x, int incx)
{
	return cblas_icamin(n,x,incx);
}

template<> ptrdiff_t warp::blas::amin<double>(size_t n, const std::complex<double> *x, int incx)
{
	return cblas_izamin(n,x,incx);
}
#endif

template<> void warp::blas::gemv<float>( size_t m, size_t n, const float* a, int lda, const float* x, int incx, float* y, int incy, const float& alpha, const float& beta, trans_e trans)
{
	cblas_sgemv(CblasRowMajor, to_blas_trans(trans), m, n, alpha, a, lda, x, incx, beta, y, incy);
}

template<> void warp::blas::gemv<double>( size_t m, size_t n, const double* a, int lda, const double* x, int incx, double* y, int incy, const double& alpha, const double& beta, trans_e trans)
{
	cblas_dgemv(CblasRowMajor, to_blas_trans(trans), m, n, alpha, a, lda, x, incx, beta, y, incy);
}

template<> void warp::blas::gemv< std::complex<float> >( size_t m, size_t n, const std::complex<float>* a, int lda, const std::complex<float>* x, int incx, std::complex<float>* y, int incy, const std::complex<float>& alpha, const std::complex<float>& beta, trans_e trans)
{
	cblas_cgemv(CblasRowMajor, to_blas_trans(trans), m, n, &alpha, a, lda, x, incx, &beta, y, incy);
}

template<> void warp::blas::gemv< std::complex<double> >( size_t m, size_t n, const std::complex<double>* a, int lda, const std::complex<double>* x, int incx, std::complex<double>* y, int incy, const std::complex<double>& alpha, const std::complex<double>& beta, trans_e trans)
{
	cblas_zgemv(CblasRowMajor, to_blas_trans(trans), m, n, &alpha, a, lda, x, incx, &beta, y, incy);
}

template<> void warp::blas::ger<float>(size_t m, size_t n, const float* x, int incx, const float* y, int incy, float* a, int lda, const float& alpha)
{
	cblas_sger(CblasRowMajor, m, n, alpha, x, incx, y, incy,a, lda);
}

template<> void warp::blas::ger<double>(size_t m, size_t n, const double* x, int incx, const double* y, int incy, double* a, int lda, const double& alpha)
{
	cblas_dger(CblasRowMajor, m, n, alpha, x, incx, y, incy,a, lda);
}

template<> void warp::blas::ger< std::complex<float> >(size_t m, size_t n, const std::complex<float>* x, int incx, const std::complex<float>* y, int incy, std::complex<float>* a, int lda, const std::complex<float>& alpha)
{
	cblas_cgeru(CblasRowMajor, m, n, &alpha, x, incx, y, incy,a, lda);
}

template<> void warp::blas::ger< std::complex<double> >(size_t m, size_t n, const std::complex<double>* x, int incx, const std::complex<double>* y, int incy, std::complex<double>* a, int lda, const std::complex<double>& alpha)
{
	cblas_zgeru(CblasRowMajor, m, n, &alpha, x, incx, y, incy,a, lda);
}

template<> void warp::blas::gerc<float>(size_t m, size_t n, const std::complex<float>* x, int incx, const std::complex<float>* y, int incy, std::complex<float>* a, int lda, const std::complex<float>& alpha)
{
	cblas_cgerc(CblasRowMajor, m, n, &alpha, x, incx, y, incy,a, lda);
}

template<> void warp::blas::gerc<double>(size_t m, size_t n, const std::complex<double>* x, int incx, const std::complex<double>* y, int incy, std::complex<double>* a, int lda, const std::complex<double>& alpha)
{
	cblas_zgerc(CblasRowMajor, m, n, &alpha, x, incx, y, incy,a, lda);
}

template<> void warp::blas::symv<float>(size_t n, const float* a, int lda, const float* x, int incx, float* y, int incy, const float& alpha, const float& beta, uplo_e upper)
{
	cblas_ssymv(CblasRowMajor, to_blas_uplo(upper), n, alpha, a, lda, x, incx, beta, y, incy);
}

template<> void warp::blas::symv<double>(size_t n, const double* a, int lda, const double* x, int incx, double* y, int incy, const double& alpha, const double& beta, uplo_e upper)
{
	cblas_dsymv(CblasRowMajor, to_blas_uplo(upper), n, alpha, a, lda, x, incx, beta, y, incy);
}

template<> void warp::blas::syr<float>(size_t n, const float *x, int incx, float *a, int lda, const float& alpha, uplo_e upper)
{
	cblas_ssyr(CblasRowMajor, to_blas_uplo(upper), n, alpha, x, incx, a, lda);
}

template<> void warp::blas::syr<double>(size_t n, const double *x, int incx, double *a, int lda, const double& alpha, uplo_e upper)
{
	cblas_dsyr(CblasRowMajor, to_blas_uplo(upper), n, alpha, x, incx, a, lda);
}

template<> void warp::blas::syr2<float>(size_t n, const float *x, int incx, const float *y, int incy, float *a, int lda, const float& alpha, uplo_e upper)
{
	cblas_ssyr2(CblasRowMajor, to_blas_uplo(upper), n, alpha, x, incx, y, incy, a, lda);
}

template<> void warp::blas::syr2<double>(size_t n, const double *x, int incx, const double *y, int incy, double *a, int lda, const double& alpha, uplo_e upper)
{
	cblas_dsyr2(CblasRowMajor, to_blas_uplo(upper), n, alpha, x, incx, y, incy, a, lda);
}

template<> void warp::blas::gemm<float>(size_t m, size_t n, size_t k, const float *a, int lda, const float *b, int ldb, float *c, int ldc, const float& alpha, const float& beta, conjtrans_e transa, conjtrans_e transb)
{
	cblas_sgemm(CblasRowMajor, to_blas_trans(transa), to_blas_trans(transb), m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);
}

template<> void warp::blas::gemm<double>(size_t m, size_t n, size_t k, const double *a, int lda, const double *b, int ldb, double *c, int ldc, const double& alpha, const double& beta, conjtrans_e transa, conjtrans_e transb)
{
	cblas_dgemm(CblasRowMajor, to_blas_trans(transa), to_blas_trans(transb), m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);
}

template<> void warp::blas::gemm< std::complex<float> >(size_t m, size_t n, size_t k, const std::complex<float> *a, int lda, const std::complex<float> *b, int ldb, std::complex<float> *c, int ldc, const std::complex<float>& alpha, const std::complex<float>& beta, conjtrans_e transa, conjtrans_e transb)
{
	cblas_cgemm(CblasRowMajor, to_blas_trans(transa), to_blas_trans(transb), m, n, k, &alpha, a, lda, b, ldb, &beta, c, ldc);
}

template<> void warp::blas::gemm< std::complex<double> >(size_t m, size_t n, size_t k, const std::complex<double> *a, int lda, const std::complex<double> *b, int ldb, std::complex<double> *c, int ldc, const std::complex<double>& alpha, const std::complex<double>& beta, conjtrans_e transa, conjtrans_e transb)
{
	cblas_zgemm(CblasRowMajor, to_blas_trans(transa), to_blas_trans(transb), m, n, k, &alpha, a, lda, b, ldb, &beta, c, ldc);
}

template<> void warp::blas::syrk<float>(size_t n, size_t k, const float *a, int lda, float *c, int ldc, const float& alpha, const float& beta, trans_e trans, uplo_e upper)
{
	cblas_ssyrk(CblasRowMajor, to_blas_uplo(upper), to_blas_trans(trans), n, k, alpha, a, lda, beta, c, ldc);
}

template<> void warp::blas::syrk<double>(size_t n, size_t k, const double *a, int lda, double *c, int ldc, const double& alpha, const double& beta, trans_e trans, uplo_e upper)
{
	cblas_dsyrk(CblasRowMajor, to_blas_uplo(upper), to_blas_trans(trans), n, k, alpha, a, lda, beta, c, ldc);
}

template<> void warp::blas::syrk< std::complex<float> >(size_t n, size_t k, const std::complex<float> *a, int lda, std::complex<float> *c, int ldc, const std::complex<float>& alpha, const std::complex<float>& beta, trans_e trans, uplo_e upper)
{
	cblas_csyrk(CblasRowMajor, to_blas_uplo(upper), to_blas_trans(trans), n, k, &alpha, a, lda, &beta, c, ldc);
}

template<> void warp::blas::syrk< std::complex<double> >(size_t n, size_t k, const std::complex<double> *a, int lda, std::complex<double> *c, int ldc, const std::complex<double>& alpha, const std::complex<double>& beta, trans_e trans, uplo_e upper)
{
	cblas_zsyrk(CblasRowMajor, to_blas_uplo(upper), to_blas_trans(trans), n, k, &alpha, a, lda, &beta, c, ldc);
}

template<> void warp::blas::syr2k<float>(size_t n, size_t k, const float *a, int lda,const float *b, int ldb, float *c, int ldc, const float& alpha, const float& beta, trans_e trans, uplo_e upper)
{
	cblas_ssyr2k(CblasRowMajor, to_blas_uplo(upper), to_blas_trans(trans), n, k, alpha, a, lda, b, ldb, beta, c, ldc);
}

template<> void warp::blas::syr2k<double>(size_t n, size_t k, const double *a, int lda,const double *b, int ldb, double *c, int ldc, const double& alpha, const double& beta, trans_e trans, uplo_e upper)
{
	cblas_dsyr2k(CblasRowMajor, to_blas_uplo(upper), to_blas_trans(trans), n, k, alpha, a, lda, b, ldb, beta, c, ldc);
}

template<> void warp::blas::syr2k< std::complex<float> >(size_t n, size_t k, const std::complex<float> *a, int lda,const std::complex<float> *b, int ldb, std::complex<float> *c, int ldc, const std::complex<float>& alpha, const std::complex<float>& beta, trans_e trans, uplo_e upper)
{
	cblas_csyr2k(CblasRowMajor, to_blas_uplo(upper), to_blas_trans(trans), n, k, &alpha, a, lda, b, ldb, &beta, c, ldc);
}

template<> void warp::blas::syr2k< std::complex<double> >(size_t n, size_t k, const std::complex<double> *a, int lda,const std::complex<double> *b, int ldb, std::complex<double> *c, int ldc, const std::complex<double>& alpha, const std::complex<double>& beta, trans_e trans, uplo_e upper)
{
	cblas_zsyr2k(CblasRowMajor, to_blas_uplo(upper), to_blas_trans(trans), n, k, &alpha, a, lda, b, ldb, &beta, c, ldc);
}

#endif
